#pragma once

namespace alt{

    //get solution
    template<size_t N,typename ValueType>//declaration
	typename variant_impl<ValueType>::solution_container get_solutions(const algebraic_equation<ValueType,N>& eq)
    {
    }

    template<typename ValueType>//linear equation
    typename variant_impl<ValueType>::solution_container
        get_solutions(const algebraic_equation<ValueType,1>& eq)
    {
		typedef ValueType value_type;
        typedef algebraic_equation<value_type,1> EqType;

        typename variant_impl<value_type>::solution_container res(EqType::order);

        res[0] = -eq[0]/eq[1];

        return res;
    }

    template<typename ValueType>//quadratic equation
    typename variant_impl<ValueType>::solution_container
        get_solutions(const algebraic_equation<ValueType,2>& eq)
    {
		typedef ValueType value_type;
        typedef algebraic_equation<value_type,2> EqType;
        typedef typename EqType::solution_type solution_type;

        const value_type& a = eq[2];
        const value_type& b = eq[1];
        const value_type& c = eq[0];

        const solution_type D = b*b - 4.*a*c;
        const solution_type sqrt_D = sqrt(D);

        typename variant_impl<value_type>::solution_container res(EqType::order);

        res[0] = (-b+sqrt_D)/(2.*a);
        res[1] = (-b-sqrt_D)/(2.*a);

        return res;
    }

    //cubic equation
    template<typename ValueType>
    typename variant_impl<ValueType>::solution_container
        get_solutions(const algebraic_equation<ValueType,3>& eq)
    {
        typedef algebraic_equation<ValueType,3> EqType;
        typedef typename EqType::value_type value_type;
        typedef typename EqType::solution_type solution_type;

        const value_type& a = eq[3];

        const value_type A2 = eq[2]/a;
        const value_type A1 = eq[1]/a;
        const value_type A0 = eq[0]/a;

        const solution_type omega(-1.0/2,sqrt(3.0)/2);

        const value_type b = -(-2.*A2*A2*A2 + 9.*A1*A2 - 27.*A0);
        const value_type _c = (A2*A2 - 3.*A1);
        const value_type c = _c*_c*_c;

        const value_type D = b*b - 4.*c;

        solution_type s0,s1,s2;

        if(D >= 0){
            const value_type sqrt_D = sqrt(value_type(D));

            const value_type s1_3 = (-b+sqrt_D)/2.;
            const value_type s2_3 = (-b-sqrt_D)/2.;

            const value_type s_s1 = s1_3>0?1:-1;
            const value_type s_s2 = s2_3>0?1:-1;

            s0 = -A2;
            s1 = s_s1*pow(std::abs(s1_3),1.0/3);
            s2 = s_s2*pow(std::abs(s2_3),1.0/3);
        }
        else{
            const solution_type sqrt_D = sqrt(solution_type(D));

            const solution_type s1_3 = (-b+sqrt_D)/2.;
            const solution_type s2_3 = (-b-sqrt_D)/2.;

            s0 = -A2;
            s1 = pow(s1_3,1.0/3);
            s2 = pow(s2_3,1.0/3);
        }

        typename variant_impl<ValueType>::solution_container res(EqType::order);      
        res[0] = (s0 + s1 + s2)/3.;
        res[1] = (s0 + omega*omega*s1 + omega*s2)/value_type(3);
        res[2] = (s0 + omega*s1 + omega*omega*s2)/value_type(3);

        return res;
    }

    ////quartic equation
    //variant_impl<complex_t>::solution_container
    //    get_solution(const algebraic_equation<complex_t,4>& eq)
    //{
    //}
}
